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We  derives  inflow  and  outflow  boundary  conditions  for  the  incompressible  Navier- 
Stokes  equations  in  cylindrical  geometries.  The  purpose  of  these  boundary  conditions  is 
to  allow  computations  in  a  finite  domain,  that  model  flow  in  an  unbounded  domain,  in  a 
way  that  the  accuracy  of  the  finite  difference  solution  is  retained,  making  the  computation 
more  efficient.  We  use  an  approach  similar  to  [9]  to  represent  the  solution  asymptotically, 
far  downstream  and  upstream,  as  a  series  expansion  which  involves  eigenvalues  and  eigen¬ 
functions.  These  eigensolutions  satisfy  certain  systems  of  ordinary  differential  equations. 
The  boundary  conditions  are  represented  by  a  family  of  differential  operators  in  a  way 
similar  to  what  was  done fin  (6]  by  Bayliss,  Gunzberger  and  Turkel.  To  demonstrate  the 
effectiveness  of  these  boundary  conditions  we  applied  them  in  numerical  computations  of 
the  incompressible  Navier-Stokes  equations  in  a  channel  with  a  step  and  in  a  pipe  with  a 
sudden  enlargement  of  the  cross  section.  To  numerically  solve  the  Navier-Stokes  equations 
we  used  a  second  order  accurate  finite  difference  schem^4^t-r^wer^ai^®i)-:,a^so  l^e  bound¬ 
ary  operators  were  approximated  using  second  order  accurate  finite  difference  formulas. 
The  numerical  results  show  the  effectiveness  and  the  increase  accuracy  obtained  by  using 
the  higher-order  boundary  conditions. 
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COMPUTATIONAL  BOUNDARY  CONDITIONS  FOR  THE  INCOMPRESSIBLE 
NAVIER-STOKES  EQUATIONS  IN  CHANNELS  AND  PIPES 

Gerardo  A.  Ache' 

1.  Introduction 

The  development  of  boundary  conditions  for  fluid  flow  in  computational  domains  is 
an  important  subject  in  engineering  applications,  e.g.  Engquist  and  Majda  [12],  Hedstrom 
[18],  Bayliss  and  Turkel  [6],  Rudy  and  Strikwerda  [29],  Han  and  Innis  [17].  In  particular 
the  determination  of  inflow  and  outflow  boundary  conditions  for  channel  and  cylindrical 
geometries  is  a  very  interesting  problem  in  computational  fluid  dynamics.  These  boundary 
conditions  are  important  since  chosen  properly  they  allow  the  size  of  the  computational 
domain  to  be  reduced  and  therefore  cut  down  the  time  of  computation.  They  also  have  a 
significant  impact  on  the  overall  accuracy  of  the  solution.  In  this  paper  we  develop  a  class 
of  inflow  and  outflow  boundary  conditions  for  the  incompressible  stationary  Navier-Stokes 
equations.  We  consider  that  class  of  infinite  domains  which  are  combinations  of  cylinders 
or  strips,  (e.g.  a  channel  with  a  step  or  a  pipe  with  a  sudden  enlargement  of  the  cross 
section),  and  for  which  conditions  for  flow  into  and  out  of  the  domain  are  needed  in  order 
to  reduce  the  finite  domain  to  a  finite  one. 
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The  unsteady  flow  of  a  viscous  incompressible  fluid  satisfies  the  Navier-Stokes  equa¬ 
tions 

~  +  (u  •  V)u  +  Vp  -  uV2u  =  f  ,  (1.1a) 

ot 

V  •  u  =  0  ,  (1.16) 

where  u(f, t)  is  the  velocity,  p(x,t)  is  the  pressure,  f(x,  t)  represents  forces  on  the  fluid, 
u  is  the  kinematic  viscosity,  x  is  in  the  domain  fi  in  RN  (N  =  2  or  3)  and  t  >  0  .  We 
here  assume  that  the  force  f  is  derivable  from  a  scalar  potential  P,  i.e.  f  =  -VP  .  Then 
we  may  rewrite  the  pressure  as  p  -f  P  .  In  addition  to  (1.1a),  (1.1b)  we  may  consider  the 
boundary  conditions 

u  =  g  on  dfl  ,  (1.1c) 

where  g  is  required  to  satisfy 

/  g-n  =  0  ,  (l.ld) 

Jan 

and  the  initial  condition 

u(f,0)  =  h(x)  .  (1.1c) 

The  Navier-Stokes  equations  are  important  since  they  describe  fluids  such  as  air  at  low 
speeds,  water,  oil,  etc.  The  mathematical  aspects  such  as  existence,  uniqueness,  and 
regularity  of  solutions  related  to  these  equations  can  be  found  in  [24]  and  [32]. 

For  the  steady  state  Navier-Stokes  equations  in  an  infinite  cylinder  with  cross  section 
C  there  is  a  relatively  simple  solution  called  Poiseuille  flow.  If  xj  measures  the  distance 
along  the  axis,  for  this  solution  the  velocity  has  the  form  ucc(x)  =  (u(t2, x3), 0. 0)  ,  and 


2 


(1.2a) 


the  pressure  p°°  =  —  Px\  ,  with  P  a  constant.  Then  u  is  solution  of  the  equation 

/  d2u  d2u  \  _  - 

U\dx\^dxl)- 

and  for  which  the  flow  condition 


J  udx^dxz  —  Q  >  0  , 


(1.26) 


is  satisfied. 

In  [l],  [4],  [20]  and  [21],  are  proven  estimates  which  show  that  the  difference  u(x)  — 
u°°(x)  decays  exponentially  in  the  axial  direction  of  a  semi-infinite  cylinder.  This  result 
is  used  to  formulate  our  boundary  conditions. 

Since  cylindrical  geometries  lead  to  Poiseuille  flow  at  infinite  distance,  (inflow  or 
outflow),  one  may  impose  at  inflow  or  outflow  the  Poiseuille  profile  as  a  computational 
boundary  condition.  This  is  often  done  in  computations  involving  the  incompressible 
Navier-Stokes  equations  in  channel  and  pipes,  the  resulting  boundary  conditions  are  of 
Dirichlet  type.  These  conditions  have  been  used  for  Navier-Stokes  calculations  and  with 
different  formulations,  e.g.  stream  function-vorticity  formulations  or  primitive  variable  for¬ 
mulations  (e.g.  [5],  [11],  [22],  [25],  [27]).  Other  type  of  computational  boundary  conditions 
are  the  specification  of  derivatives  of  the  velocity  or  the  stream  function  and  vorticity  at 
inflow  or  at  outflow,  (see  [19],  [23],  [26]).  For  example,  in  [16]  Greenspan  performed  Navier- 
Stokes  calculations  using  Poiseuille  flow  as  an  inflow  boundary  condition,  but  at  outflow  he 
combined  the  vorticity  and  the  stream  function  in  a  way  that  makes  the  pressure  constant 
at  the  outflow  boundary.  He  also  carried  out  computations  using  Poiseuille  flow  at  both 
boundary  positions  and  he  pointed  out  that  both  formulations  were  numerically  equiva¬ 
lent.  In  [13]  the  effect  of  some  downstream  boundary  conditions  in  a  channel  with  a  step 


are  discussed  for  both  the  stationary  and  time  dependent  incompressible  Navier-Stokes 
equations. 

Since  the  Poiseuille  profile  occurs  only  in  the  limit  at  infinite  distances,  in  order  to  con¬ 
struct  better  boundary  conditions  we  may  consider  the  asymptotic  nature  of  the  solution 
in  the  manner  of  Bramley  and  Dennis  [9].  We  assume  that  the  solution  can  be  regarded  at 
large  distances  as  Poiseuille  flow  plus  a  perturbation  which  decays  exponentially  fast.  This 
perturbation  is  then  expanded  in  terms  of  an  eigenfunction  series  involving  eigenvalues 
and  eigenfunctions  which  satisfy  certain  system  of  ordinary  differential  equations  [l],  [8], 
[9],  [10],  [15]  and  [33].  Using  this  approach  we  may  construct  better  boundary  conditions 
by  introducing  an  operator  which  annihilates  the  first  few  terms  in  the  expansion.  This 
method  is  closely  related  to  that  presented  by  Bayliss,  Gunzburger,  and  Turkel  [5]  for  the 
Helmholtz  equation,  where  the  concept  of  asymptotic  expansion  was  applied  to  construct 
far-field  boundary  conditions,  see  also  [7],  Similar  to  [5],  the  operators  which  define  our 
boundary  conditions  are  linear  differential  operators  and  they  may  involve  derivatives  of 
order  greater  than  the  order  of  the  differential  equations.  These  boundary  conditions  are 
different  from  the  usual  conditions  used  in  incompressible  Navier-Stokes  calculations  which 
involve  only  the  variables  and  first  derivatives.  For  linear  systems  of  elliptic  equations  it 
is  known  (see  [3])  that  it  is  possible  to  prescribe  boundary  conditions  defined  in  terms  of 
linear  differential  operators  involving  derivatives  higher  than  the  order  of  the  system  and 
still  have  the  problem  be  well-posed. 


In  [9]  Bramley  and  Dennis  formulated  inflow  and  outflow  boundary  conditions  in 
channel-like  geometries  which  are  similar  to  our  first-order  boundary  conditions,  but  ap- 
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plied  to  the  stream  function  formulation  of  the  Navier-Stokes  equations.  They  suggest 


that  their  conditions  may  be  suitable  for  numerical  work  involving  the  incompressible 


Navier-Stokes  equations. 


To  implement  the  boundary  conditions  we  used  finite  difference  formulas  for  the 


derivatives  as  efficiently  as  possible  (these  approximations  are  of  the  same  order  of  ac¬ 


curacy  as  the  order  of  the  numerical  scheme  which  is  second  order  accurate  (Strikwerda 


[30])).  The  numerical  results,  as  discussed  in  section  4,  show  the  effectiveness  and  the 


increase  of  accuracy  obtained  by  using  the  higher  order  boundary  conditions. 


2.  Development  of  Inflow  and  Outflow  Boundary  Conditions 


Similar  to  Bramley  and  Dennis  [9],  we  assume  that  the  solution  to  the  Navier-Stokes 


equations  (u(x),p(x))  can  be  represented  asymptotically,  far  upstream  or  far  downstream, 


in  the  form 


u(x)  =  u°°(y)  +  Wn(y)exp(-ABxi) 


where  x  =  (i1?Z2 ,  ••*,xjv)  =  (xj,y)  ,  and  u°°(y)  is  the  Poiseuille  profile  associated  with 


the  infinite  cylinder.  The  set  j  Wn(y)  exp(-A„xi)  j  is  assumed  to  be  a  complete  set 
of  eigenfunctions  with  |wn(y),An|  being  solution  to  an  eigenvalue  problem  which  is 


described  in  §3.  The  eigenvalues  An  satisfy, 


An  |<|  An4-i  |  for  n  =  1,2, 


Eigenvalues  with  negative  real  part  are  associated  to  the  asymptotic  solution  upstream 


(i.e.  ij  negative)  and  those  with  positive  real  part  with  the  solution  downstream  (i.e.  xj 


positive). 
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To  construct  an  artificial  boundary  at  x\  =  d  (inflow  or  outflow)  we  observe  that  from 
(2.1)  we  get  the  following  relation, 

3 


dx i 


(u  -  u°°)  +  Ai(u  -  u°°)  =  0[exp(-A2d)j 


(2.3) 


t  I  *,=<* 

i.e.  the  operator  djdx\  4-  Ai  annihilates  the  first  term  in  the  expansion  (2.1).  So  we  may 
define  the  following  boundary  condition  at  outflow  or  inflow,  for  d  large, 

d 


(u  -  u°°)  +  Ax(u  -  u°°)  =0 

dx\  ii1=d 


(2.4) 


A  more  accurate  boundary  condition  is  obtained  by  annihilating  the  first  two  terms  in 
(2.1),  and  has  an  error  of  magnitude  0[exp(-A3d)l  .  The  operator  is 


<£  +  a‘»£  +  a’>  • 


(2.5) 


In  case  that  Aj  is  complex  with  A2  =  ^  (2.5)  becomes 


— +  2fie(Ai)^-+|A1 


(2.6) 


To  develop  even  more  accurate  conditions  at  the  boundary,  we  consider  the  sequence  of 


operators 


-  ri(aI1  +a*)  • 


»=i 


(2.7) 


A  straightforward  calculation  gives  that  £m(v)  =  0  ,  for  any  vector  function  v  of  the  form 


v(£)  =  ^f(y)exp(-Akxj) 


(2.8) 


k  —  1 


Using  the  representation  for  u  in  (2.1)  we  obtain  that 


Bm{  u-u°°)  =  0[e\p{-Xm^idy  . 

ii  =d 


(2.9) 


That  is  for  any  integer  m  >  0  the  differential  operator  Bm  annihilates  the  first  m  terms  in 
the  asymptotic  expansion  (2.1).  Therefore  we  can  use  this  family  of  operators  as  a  way  of 
matching  asymptotically  the  solution  up  to  the  first  m  terms.  We  may  therefore  impose, 
at  inflow  or  outflow,  the  following  boundary  conditions 

Bm{  u-u°°);ii=d  =  0  ,  (2.10) 

which  is  accurate  at  the  boundary  with  an  error  of  order  O  |exp(-Am+1d)  ]  .  We  may 
define  a  zero-order  boundary  condition  replacing  in  (2.10)  Bm  by  the  identity  operator. 
Then  for  m  =  0  ,  the  boundary  condition  of  this  family  is  of  a  Dirichlet  type  which  is  the 
most  common  boundary  condition  used  at  inflow  or  outflow  in  engineering  computations. 
However,  one  can  expect  a  significant  improvement  in  the  efficiency  of  the  computations 
by  using  the  higher-order  boundary  conditions  rather  than  the  Dirichlet  condition  because 
of  faster  decrease  of  the  error  at  the  inflow  or  outflow  boundary.  The  result  will  be  a 
considerable  gain  in  computer  storage  and  running  time  as  well,  due  to  a  reduction  in 
the  size  of  the  computational  domain.  Notice  that  since  the  operators  which  define  these 
boundary  conditions  involve  higher-order  derivatives  we  also  can  expect  an  increase  of  the 
difficulties  in  the  implementation  of  these  boundary  conditions  as  the  order  increases.  As 
we  will  show,  these  difficulties  in  the  implementation  are  usually  offset  by  the  substantia! 
increase  in  accuracy  and  efficiency. 

3.  The  Eigenvalue  Problems 

The  eigenpair  ( Wn(y),  An)  are  obtained  from  an  eigenvalue  problem  that  results  when 
we  seek  asymptotic  solutions  of  the  Navier-Stokes  equations  in  the  form 


p(i)  =  p“(ii)  +  R  '{(sOexpf-Ai,)  . 


(3.16) 


For  flow  in  a  channel,  neglecting  nonlinear  terms,  these  eigenvalue  problems  can  be  written 


d2W, 
dy 2 


-fl{u(y) AH',  -  —W2)  -  \q-\2\\\  , 


=  A W,  , 


dW  y 

dy 


(3.2a) 


(3.26) 


(3.2c) 


with  boundary  conditions 


W,(±1)  =  W'2(±1)=0 


(3. 2d) 


Similarly  for  axi-symmetric  flow  in  a  pipe 


r.  W, 

2  =  \W2  -  —  , 


d2W2 


d{lix,  x.-./.m*  1  1<m"2  \  2i 


-  Au(r)W'2}  -  -  A2W2  -  A*  , 

dr-  dr  r  dr 

^  =  #u(r)AW',  -r  A^  +  A 2 IV,  , 
dr  dr 


(3.3a) 


(3.36) 


(3.3c) 


with  boundary  conditions, 


W  ,(0)  -  — — 1  (0)  =  0  and  IV, (1)  =  W2(l)  =  0  , 
dr 


(3.3  d) 


where  u  is  the  parabolic  profile  in  Poiseuille  flow,  as  defined  in  (1.2a),  (1.26)  and  R  =  ® 
is  the  Reynolds  number.  Results  concerning  the  numerical  solution  of  these  eigenvalue 
problems  can  be  found  in  [l],  [8],  [9],  {10],  [15]  and  !33],  for  flow  in  a  channel  and  in  [1], 
for  axi-symmetric  flow  in  a  pipe,  also  there  are  even  and  odd  eigenvalues  associated  to  the 
channel  flow  problem.  In  this  paper  we  compute  solutions  of  the  Navier-Stokes  equations. 
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in  a  channel  and. a  pipe,  using  inflow  and  outflow  boundary  conditions  of  the  form  (2.10), 
where  the  eigenvalues  A,  in  Bm  are  determined  from  the  eigenvalue  problems  (3.2)  and 
(3.3)  respectively. 

4.  Numerical  Implementation 

An  important  ingredient  needed  to  test  these  boundary  conditions  is  a  numerical 
scheme  to  solve  the  Navier-Stokes  equations.  The  method  we  used  here  is  based  on  a 
second-order  finite  difference  scheme  to  solve  the  Stokes  and  Navier-Stokes  equations  ;29\ 
[30].  This  method  deals  with  primitive  variables  formulation,  i.e.,  it  uses  the  velocity  and 
the  pressure  as  dependent  variables.  Then  the  laplacian  and  convection  terms,  given  in 
conservative  form,  are  discretized  using  standard  centered  difference  formulas,  while  the 
pressure  gradient  and  the  continuity  equation  are  discretized  using  regularized  centered 
difference  formulas  [29],  e.g. 


«  <W  -  ^h'16I-6l+p  , 

(4.1a) 

p-\h26y-6l+p  . 

(4.16) 

& u  1  r  1,2  r  c  2 

~  t>ioui  _  h  oI4_oi_ui  , 

ox  6 

(4.1c) 

du 2  .  1  ,  2  c  r  2 

Qy  —  U  y0  ^  2  g  ^  Vy4-Oj(_U2  , 

(4.1d) 

where  (<5xo,<5yo),  (6z-,6y_),  (6z+,6y+)  are  centered,  backward,  and  forw-ard  differences  in 
the  x  and  y  direction,  respectively.  Formulas  (4.1)  are  used  in  the  pressure  gradient  and 
the  continuity  equation.  To  determine  the  pressure  on  the  boundaries  cubic  interpolation 


formulas  are  used,  i.e. 


(4.2c) 


-V  *'■  ,*v 

POJ  =  3(Plj  -  P2,>)  —  Pz,j  i 


P.,0  -  3(ptii  -  p,t 2)  -  Pi, 3 


(4.26) 


To  solve  the  resulting  finite  difference  equations  an  extended  S.O.R.  iterative  proce¬ 
dure  [30j  was  used.  This  numerical  method  is  very  efficient  and  accurate.  It  has  been  used 
to  compute  flow  in  a  spinning  and  coning  cylinder  j 31 ; . 

To  implement  the  boundary  conditions,  we  construct  numerical  approximations  to 
these  derivatives  of  order  0(h2),  where  h  is  the  grid  size  in  the  axial  direction.  These 
approximations  do  not  destroy  the  entire  accuracy  of  the  interior  scheme.  We  provide  the 
finite  differences  formulas  used  to  approximate  the  differential  operators  for  each  boundary 
condition,  e.g. 

— First-order  : 

av_  =  CAT-C-,  +0{h2)  (4.4) 

&xl  TV  — 1/2  ^ 


— Second-order: 


a2 


u 


dx2 


un  -  2u/v_i 


u  \  —  2 


N  -  1 


h 2 


0(h 2) 


(4.5a) 


— Third-order  : 


du  I 

&x\  :a'_i 


2 h 


0(h 2) 


(4.56) 


dzu 

dx3 


N-  3.  2 


u.v  _  3(u.v-i  a.v-o)  —  w.v-3 

h3 


0{h 2) 
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(4.6a) 


d2u  UN  -  {us-\  +  UN-2)  +-UN-3  ,  ^/.2x 

dx 2  '  '  •  ~  2h2  T  ^  ’ 

axl  .V  — 3/2 


(4.66) 


UN_i  -  UJV-2  ,  2v 


I  TV — 3/2 


+  0{h2)  . 


(4.6c) 


— Fourth-order  : 


d4u  UN  -  4uN_!  -r  6uA'_2  -  4uA  _3  4-  U.v_4  .  2 


<9r4  ! 

axi  i a--: 


0(/i2)  ,  (4.7a) 


<9r3  ' 
axi  in-: 


«Af  —  2(ujv-l  -  U A’  — 3)  —  tiN-4 


0(/r2) 


(4.76) 


d2U  Ua-1  -  2ua  — 2  +  «A-3  ,  ~,,2X 

=  -  -  +0[h)  , 

axl  N- 2  n 


(4.7c) 


Ua-1  -  «A'-3 


^X1  A-5 


Hh  0(h2)  . 


(4-7-0 


5.  Test  Problems 

In  order  to  test  these  boundary  conditions  we  used  two  problems  as  test  problems. 
The  first  problem  corresponds  to  an  incompressible  viscous  flow  in  a  channel  with  a  step 
(see  Figure  1)  which  can  be  regarded  as  the  union  of  two  strips,  one  upstream  of  the  form 
Ri  =  (5,1)  x  (— oo.Oj  and  the  strip  downstream  of  the  form  R2  =  (0,1)  x  ;0,oo)  . 


*‘.v 


The  Poiseuille  solution  associated  with  the  infinite  strip  having  the  same  cross  section 


as  R j  is  given  by 


«(y)  =  48(y  -  ^)(1  -  y)  , 

(5.1a) 

p°°(x)  = -^x  + C,  , 

(5.16) 

where  C\  is  a  constant. 

Similarly  the  Poiseuille  solution  associated  with  the  infinite  strip  with  same  cross 

♦ 

section  as  Ro  is  given  by 

«(y)  =  6y(l  -  y)  ,  (5.2a) 

P°°(x)  = -^x  +  C2  .  (5.26) 

To  solve  the  Navier-Stokes  equations  we  imposed  boundary  conditions  at  inflow  and 
outflow  of  the  form  (2.8),  and  the  non-slip  condition  at  the  walls.  The  eigenvalues  in  the 
expansion  (2.1)  are  obtained  by  solving  the  eigenvalue  problems  described  in  §3. 

The  second  test  problem  consists  of  a  pipe  with  a  sudden  enlargement  on  the  cross 
section,  the  fluid  motion  is  considered  axi-symmetric.  Computations  are  made  only  in 
half  of  the  domain,  so  the  computational  domain  is  similar  to  that  for  the  channel  flow 
problem. 

The  Poiseuille  solutions  associated  which  each  part  of  the  domain  are  given  by, 

u(r)  =  4(1  -  4r2)  ,  (5.3a) 


P~(*) 


Ci  . 


(5.36) 


for  the  pipe  upstream,  i.e.  R\  -  (0,  I)  x  (-oo,0  and 


for  the  pipe  downstream,  i.e.  R2  =  (0, 1)  x  (0, 00). 

These  two  problems  have  been  the  subject  of  many  numerical  studies  in  the  last  20 
years,  e.g.  [11],  [16],  [19],  [22],  [23],  [25],  [26],  [27]. 

6.  Numerical  Results 


Figure  1:  A  Channel  with  a  step. 

We  have  done  numerical  computations  which  include  the  two  test  problems  for  flow- 
in  a  channel  and  flow  in  a  pipe.  These  two  problems  retain  the  full  non-linear  character 
of  the  Navier-Stokes  equations.  Boundary  conditions  of  order  0,  1,  2,  3  and  4  have  been 
applied  at  outflow  for  the  channel  flow-  problem  and  boundary  conditions  of  order  0,  1  and 
2  have  been  applied  at  outflow  for  the  pipe  flow  problem.  The  specification  of  Poiseuille 
flow  (zero-order  boundary  condition)  was  applied  at  inflow  for  both  problems.  The  test 


positions  where  the  computational  boundaries  were  specified  were  at  x  =  d  =  —  1  (inflow) 
and  x  —  d,  —  A  (outflow).  The  computations  have  been  performed  for  Reynolds  numbers 
R  in  the  range  0  <  R  <  50,  and  solutions  were  obtained  in  the  largest  domain  {d  =  4), 
using  the  grids  110  x  21,  140  x  29,  and  170  x  35  at  R  =  30,  40,  and  50,  respectively,  for  the 
channel  flow  problem.  Similarly  for  the  pipe  flow  calculations  the  grid  sizes  were  50  x  21, 
.70  x  21,  and  110  x  21  at  R  =  10,  20,  and  30,  respectively.  The  iterative  procedure  was 
stopped  "when  the  residuals  of  the  difference  satisfied  a  given  tolerance  of  .5  x  10-5  in  the 
largest  domain  and  .5  x  10-4  for  computations  in  the  shorter  domain  (i.e.  d  <  4).  All  the 
computations  were  done  on  the  VAX  11/780  at  the  Mathematics  Research  Center  at  the 
University  of  Wisconsin-Madison. 

The  numerical  study  included  the  following  aspects  related  to  the  boundary  conditions. 
For  a  fixed  domain,  we  considered  the  effect  of  these  boundary  conditions  as  the  Reynolds 
number  was  increased.  For  a  fixed  Reynolds  number,  we  considered  the  effect  of  these 
boundary  conditions  as  the  domain  was  reduced. 

In  Figures  2.1a  to  2.3b  we  display  graphs  corresponding  to  the  velocity  field  at  the 
center  line  of  R2,  for  the  channel  and  pipe,  downstream,  i.e,  at  the  line  T  =  {(x,  y)/y  = 
I,a  <  x  <  d},  w'here  a  is  greater  than  zero  (i.e,  away  from  the  step).  These  graphs  show 
how  the  solution  is  affected  when  boundary  conditions  of  different  orders  are  applied.  In 
the  channel  flow  problem  at  R  equal  to  30  Figures  2.1a  to  2.1b  show  the  velocity  field 
corresponding  to  the  domain  at  different  distances  downstream  (i.e,  d  is  4  and  d  is  1.5). 
these  solutions  correspond  to  boundary  conditions  of  order  0.2.3,  and  4.  Since  at  that 
Reynolds  number  the  first  and  second  eigenvalues  in  the  expansion  (2.1)  are  a  complex 


conjugate  pair  (even)  and  the  third  and  fourth  eigenvalues  are  real  (odd  and  even),  see 
Table  1,  it  is  not  appropriate  to  use  a  first-order  boundary  condition.  In  Figure  2.1a  we 
notice  that  the  solutions  for  the  different  boundary  conditions  are  close  together,  except 
near  the  outflow  boundary  where  the  solution  corresponding  to  the  zero-order  boundary 
condition  presents  the  largest  deviation  from  solutions  obtained  using  the  higherr-order 
conditions.  As  the  domain  is  reduced  from  d  of  4  to  1.5  the  solution  corresponding  to  the 
zero-order  condition  separates  from  the  solutions  corresponding  to  the  higher-order  condi¬ 
tion,  these  solutions  remaining  close  together.  We  computed  the  absolute  and  percentage 
errors  for  the  difference  between  the  solution  in  the  largest  domain  and  solutions  in  a 
shorter  domain,  for  each  boundary  condition,  these  errors  were  measured  using  the  dis¬ 
cretized  L2*norm5  i-e-i  if  Uh  ‘s  the  numerical  solution  in  the  largest  domain,  corresponding 
to  boundary  condition  of  order  fc,  and  Vk  is  the  numerical  solution  in  a  shorter  domain, 
then  the  absolute  and  relative  error  are  defined  as  Afc  :=  |jf/£  -  Vfcfc|j2  and  Ak/\\U£\\2  re¬ 
spectively,  where  ||  || 2  :=  MS/  Si  U ij )1//2>  and  the  sum  is  taken  over  all  the  grid  points. 

We  have  found  that  at  R  equal  to  30,  the  numerical  solution  obtained  using  the  zero-order 
condition  in  the  domain  with  outflow  test  positions  at  d  equal  to  1.5  present  percentage 
errors  of  about  12%,  while  for  solutions  corresponding  to  the  higher-order  condition  these 
percentage  errors  were  less  than  1%. 

At  Reynolds  number  40  and  50  the  first  two  eigenvalues  in  (2.1)  are  real  (odd  and 
even)  and  the  third  eigenvalue  is  complex  (even)  and  thus  we  applied  boundary  conditions 
of  order  0, 1,2  and  4.  We  notice  that  at  R  equal  to  40  the  gap,  near  the  outflow  boundary, 
between  the  numerical  solutions  corresponding  t,o  the  zero-order  condition  and  solutions 


corresponding  to  the  higher-order  condition  is  bigger  that  for  those  at  R  equal  to  30. 
For  d  <  3  a  numerical  solution  using  the  first-order  boundary  condition  could  not  be 
obtained,  and  for  d  equal  to  1.5  solutions  were  obtained  only  for  boundary  conditions 
of  order  2  and  4.  For  R  equal  to  50  and  d  equal  to  4,  solutions  were  obtained  for  all 
boundary  conditions.  At  d  equal  to  3  solutions  were  obtained  for  the  first,  second  and 
fourth  order  boundary  conditions,  at  d  equal  to  2  solutions  were  obtained  using  the  second 
and  fdftjth  order  boundary  condition  and  at  d  equal  to  1.5  a  solution  was  obtained  only 
for  the  fourth-order  boundary  condition.  For  these  Reynolds  numbers  the  largest  absolute 
and  percentage  errors  correspond  to  numerical  solutions  obtained  using  the  zero-order 
boundary  condition,  also  the  absolute  and  percentage  errors  for  solutions  corresponding  to 
the  higher-order  condition  become  smaller  as  the  Reynolds  number  is  increased,  that  is  due 
to  the  grid  size  being  finer.  These  grid  sizes  were  necessary  in  order  to  get  a  solution  using 
the  zero-order  boundary  condition.  However,  solutions  on  a  coarser  grid  can  be  obtained 
using  the  higher-order  boundary  conditions.  As  we  increased  the  Reynolds  number  from 
R  equal  to  50  to  R  equal  to  60  we  could  not  get  a  solution  using  the  zero  or  first  order 
condition  on  the  coarser  grid,  how'ever  as  the  grid  was  refined  a  solution  was  obtained 
using  the  first-order  condition  but  not  the  zero-order  condition. 

For  the  pipe  flow  problem  we  present  the  same  type  of  results  as  for  the  channel  flow 
problem.  In  this  case  the  effect  of  using  the  zero-order  boundary  condition  or  a  higher- 
order  condition  is  more  notable,  notice  also  that  the  eigenvalues  for  this  case  are  smaller 
in  magnitude  than  those  for  the  channel  flow  problem  (Tables  1.2).  For  this  problem  we 
applied  boundary  conditions  of  order  0  and  1  at  Reynolds  numbers  R  equal  to  10.  20  and 


boundary  conditions  of  order  0,  1  and  2  at  R  equal  to  30  ,  these  conditions  correspond  to 
real  eigenvalues  in  the  expansion  (2.1).  Similar  to  the  channel  flow  problem,  we  computed 
absolute  and  percentage  errors  for  the  different  boundary  conditions.  In  this  case  the 
percentage  errors  corresponding  to  the  zero-order  boundary  condition,  for  the  pipe  flow 
problem,  are  much  bigger  than  those  for  the  channel  flow  problem  (see  Table  3).  It  is 
then  possible  to  appreciate  the  effect  of  these  boundary  conditions  in  the  stream  function 
and  pressure  contour  plots.  In  Figures  3. a  to  3.d  we  present  stream  function  and  pressure 
contour  plots  for  different  boundary  conditions  at  different  test  positions.  Notice  that 
near  the  outflow  boundary  in  the  pressure  contour  plots  corresponding  to  the  zero-order 
condition  (Figures  3. a,  3.c)  there  are  contours  that  can  not  be  seen  for  the  pressure  contour 
plots  corresponding  to  the  higher-order  condition.  Similarly  the  stream  function  contour 
plots  associated  to  the  zero-order  condition  (Figures  3. a.  3.c)  presents  some  oscillations 
near  the  outflow  boundary  while  the  stream  function  contour  plots  corresponding  to  the 
higher-order  condition  do  not  present  those  oscillations.  These  oscillations  in  the  pressure 
and  in  the  stream  function  contour  plots  are  the  result  of  the  zero-order  boundary  condition 
being  applied.  Notice  that  in  the  stream  function  contour  plots  the  oscillations  are  not 
as  big  as  for  the  pressure  contour  plots,  that  is  because  we  computed  the  stream  function 
integrating  the  velocity,  using  the  trapezoidal  rule,  and  the  integration  acts  to  smooth  the 


7.  The  Pressure  Corner  Singularity 

Analytically,  it  can  be  shown  ([l],  [2])  that  for  plane  flow  the  pressure,  near  a  re¬ 
entrant  corner,  behave  as 

p  «  A,r-°-456  +  O(r~0091)  ,  (7.2) 

where  A\  is  a  constant,  and  the  pressure  is  given,  near  the  corner,  in  polar  coordinates, 
with  the  origin  at  the  re-entrant  corner. 

Since  the  geometries  considered  in  this  paper,  for  the  two  test  problems,  present 
re-entrant  corners  and  the  pressure  tends  to  be  unbounded  as  the  re-entrant  corner  is 
approached  we  need  to  provide,  to  the  finite  difference  scheme  finite  pressure  values  at 
such  corners,  as  efficient  as  possible,  to  evaluate  the  pressure  gradient  at  adjacent  grid 
points.  In  our  situation,  since  the  numerical  scheme  requires  values  of  the  pressure  at  the 
solid  wall  and  since  these  values  are  computed  using  extrapolation  formulas  (as  described 
in  §4  ),  we  may  compute  the  pressure  at  the  corner  in  two  different  ways  using  such 
extrapolation  formulas  for  each  of  the  two  directions.  Assuming  that  the  re-entrant  corner 
is  located  at  (xo,t/o)>  we  have 

Pc,x  ■=  p{xo,yo) 

=  3 (p(x0  +  hx,y0)  -  p{x0  +  2 hx,y0))  +  p{x0  +  3 hz,y0)  ,  (7.1a) 

and 

Pc,y  ■=  p(*o,yo) 

=  3(p(x0,y0  +  hy)  -  p(x0,yc,  +  2hy))  p(x0,y0  +  3/iy)  .  7.16) 

Here  hz,hy  are  the  grid  sizes  in  the  x  and  y  directions  respectively.  Notice  that  since  these 
values  of  PCtX  and  Pc  v  do  not  need  to  be  equal  it  suggests  that  the  pressure  at  a  re-entrant 


corner  may  be  regarded  as  a  double  valued  function.  Therefore  to  evaluate  the  pressure 
gradient  at  adjacent  grid  points  we  use  Pc<1  to  compute  the  approximation  to  dpjdx  at 
(xo  +  hx,y o)  and  Pc>y  to  compute  the  approximation  to  dp/dy  at  (xo,yo  +  hy). 

The  pressure  strategy  described  above  works  well  in  computation  at  small  and  mod¬ 
erate  Reynolds  numbers  and  accurately  models  the  pressure  singularity,  as  it  is  shown  in 
[1]  and  [2]. 


8.  Conclusion 

We  have  developed  a  family  of  local  boundary  conditions  for  the  incompressible 
Navier-Stokes  equations  in  channel  and  pipes.  The  boundary  conditions  are  constructed 
using  the  fact  that  the  solution  possesses  an  asymptotic  expansion  of  the  form  (2.1).  Nu¬ 
merical  results  demonstrate  the  effectiveness  of  these  boundary  conditions.  These  results 
are  summarized  as  follows: 

1-  For  a  given  domain,  a  given  grid  size,  and  a  given  boundary  condition,  as  the  Reynolds 
number  is  increased  a  solution  may  or  may  not  be  obtained.  In  case  the  solution  is  not 
obtained  either  the  finite  difference  solution  does  not  exist  or  the  iterative  method  fails  to 
compute  the  solution.  As  the  mesh  is  refined  a  solution  may  be  obtained. 

2-  For  a  given  domain,  a  given  grid  size,  and  a  given  Reynolds  number,  as  the  order  of  the 
boundary  condition  increases  the  accuracy  of  the  finite  difference  solution  improves. 

3-  For  a  given  grid  size,  a  given  Reynolds  number,  and  a  given  boundary  condition,  as  the 
domain  size  increases  the  accuracy  increases. 

4-  For  a  given  domain,  a  given  Reynolds  number,  and  a  given  boundary  condition,  as  the 
grid  size  decreases  the  accuracy  increases. 
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These  results  demonstrate  that  the  use  of  the  higher-order  boundary  conditions  can 
improve  both  the  efficiency  and  the  accuracy  of  fluid  flow  calculations. 
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TABLE  1 


Eigenvalues  from  the  expansion  (2.1), for  the  channel  flow  problem 


“ 

R 

inflow 

Ai 

outflow 

Ai 

outflow 

a2 

outflow  J 

A3  I 

30 

-6.33792 

1.9842— il. 2012 

2.48996 

3.372034 

40 

-5.95994 

1.86981 

1.89935 

2.1503— il. 1312 

50 

-5.69344 

1.31807 

1.49834 

2.0832— il. 2186 

60 

-5.49035 

0.95132 

1.25008 

1.9780— il  .2161 

TABLE  2 


Eigenvalues  from  the  expansion  (2.1),  for  the  pipe  flow  problem 


R 

inflow 

Ai 

outflow 

Ai 

outflow  | 

a2 

10 

-11.2963-i2.1544 

2.84871 

4.2920-fil. 42646 

20  . 

-9.89212 

1.23859 

3.9266— iO. 7599  ! 

30 

-9.38262 

1.02639 

2.83212 

TABLE  3 


Absolute  and  percentage  errors  betwwen  the  solution  in  the  largest 


domain  and  solutions  in  a  shorter  domain  (pipe  flow) 


R 

— — _l 

d 

' 

B.C 

% 

AT*) 

%  | 

10 

3 

0 

0.00048 

1.28 

0.00058 

0.23  1 

2 

0.00566 

13.10 

0.00594 

2.72 

1.5 

0.01779 

38.72 

0.01697 

8.89 

3.0 

1 

0.00028 

0.75 

0.00059 

0.23 

2.0 

0.00081 

1.87 

0.00209 

0.95 

1.5 

0.00190 

4.13 

0.00438 

2.29 

20 

3.0 

0 

.  0.01095 

36.06 

0.00869 

3.92 

2.0 

0.03650 

115.63 

0.02506 

13.33 

3.0 

1 

0.00080 

2.64 

0.00108 

0.49 

30 

3.0 

0 

0.03200 

128.68 

0.01989 

9.83 

3.0 

2 

0.00001 

0.05 

0.00003 

0.02 

The  radial  velocity  component 


*  The  axial  velocity  component 
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FIGURE  2.1b  :  Channel  flow  problem.  The  transverse  velocity  component 


centerline  for  R  =  10  Boundary 


0.2461 


FIGURE  2.2c  :  Pipe  How  problem.  The  radial  velocity  component  on 
centerline  for  R  =  30  Boundary  conditions  of  order  0,1  and  2  (o,  x,d). 


